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ABSTRACT 
Let Xyoeeeaky have a joint k-variate normal distribution with zero 


means, common unknown variance of and known correlation matrix (o55)> where 


O57 op for all i#j. Let 5° 


that vse/oe has a chi-squared distribution with v degrees of freedom. New 


be distributed independently of the x. such 


tables with wider coverage and more accuracy that the previously published 
X. 
ones are given for the percentage points of Y = max . Some basic 
l<ix<k 
theoretical results are given in Section 2. The next section describes 
Hartley's results for approximating the distribution function of Y. Besides 
a brief review of existing tables (Section 4), the paper discusses the 


construction of new tables based on Hartley's results (Section 5) and some 


specific applications (Section 6). 
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1. INTRODUCTION 


Let X .,X, have a joint k-variate normal distribution with zero 


De rk 
means, a common unknown variance a and a known correlation matrix 

Q = (p55): Also, let 5¢ be a random variable independent of the X,'s 
such that vse/oe has a chi-square distribution with v degrees of freedom. 
Then the joint distribution of ty otos---stys where t; = X./s, is a 
(central) multivariate t distribution which has the density (see Cornish 


{2] or Dunnett and Sobel [4]) 
vtk 


i -(-=>-) 
(1.1) (tyre eeaty) = ne re 01 oe) Jay stity] eo 


vt 


7 = ar! 
where A = (a; 5) Q. 


We are interested in the distribution of 


max X; 
(1.2) Y = max t; = r 


Of course, the percentage points of the distribution of Y are the same as 
the corresponding equi-percentage points of the k-variate t distribution 
in (1.1). These percentage points are needed in several problems of 

statistical inference such as tests for the multiple comparisons of means 


(or mean vectors) against one-sided alternatives under a general analysis 


*This research was supported by the Office of Naval Research Contract 
N00014-75-C-0455 at Purdue University. Reproduction in whole or in part 
is permitted for any purpose of the United States Government. 


of variance (or multivariate analysis of variance) model, simultaneous 
confidence intervals for multiple comparisons between several treatment 
means and a control mean, and selection and ranking procedures for selecting 
the best population and for selecting populations that are better than a 
control. So it is not surprising that several authors have provided tables 
of the percentage points of the distribution of Y. These tables differ 
with respect to the ranges of the various parameters, the accuracy of the 
values, the probability levels chosen, and the numerical methods employed 
in constructing these tables. 

Besides providing new tables of these percentage points with wider 
coverage and more accuracy, our emphasis in this paper is on Hartley's 
method of approximating pr{Y < y} by a series in which the leading term 
corresponds to the case of v = ~ (o° known) and the successive terms are 
in powers of oe Hartley [11] obtained this as the solution of a certain 
differential-difference equation. His method is found very useful; however, 
no published tables of the percentage points of Y until] now were based on 
this method. 

In this paper, we discuss some basic theoretical results (Section 2) 
and Hartley's results with an additional correction term (Section 3). Besides 
reviewing briefly the existing tables (Section 4), we discuss the construction 


of the new tables (Section 5), and some specific applications (Section 6). 


2. BASIC THEORETICAL RESULTS 
In considering the distribution of Y we assume without loss of 
generality that o = 1. Further, we only consider the case where P ij =o > 0 
for all i # j. In this case, it is well-known (Dunnett and Sobel [4]) that 


the random variables Ky oXooeee oXy can be generated from k+1 independent standard 


normal variables Zoreqeees ody by setting 
Beatie 3 
(2.1) uF = (1-p) z; - / Zo» i= 1,...,k. 


Using this representation, we can write the distribution function 
FL (ws co) of W = max X in the form 
” owe t 
a + 
(2.2) FL (ws 9) = f CEP) o(t)at 
k 
-2 (1-9)? 
where »(-) and ¢(-) denote the density and the cdf of a standard normal 
variate. By a similar argument, the distribution function Gy, (ys 0) 


of Y is given by 
(2.3) G (ys e) [Jee HQC, (u)dt du 
0 - 


where h (+) is the density function of s. In the sequel, unless the 
context demands us to be more specific, we will simply write G (y) and 
F(y) instead of Gs o) and F.(y; o), respectively. It is a well-known 
fact that lim G (y) = F(y). 

yo 


The moments of Y can be written as a product of two expectations. 


Let ZK] = max(Z1,...sZ)). Forr<  v, 


(2.4) ELV] = Els ]e{(vT-5 Zy4 - vb Zq)"] 
real 


-Vr 


[r/2] 
ni 03) 3 ae (1 0)? 
jz 


ake r-2j° 


where [r/2] is the largest integer less than or equal to r/2, ae i is the 
% 


ith moment of the largest of k independent standard normal variables, and 


(vis 
(2.5) = E[s*] - —4 : 
: (3°72) 


3. HARTLEY'S METHOD 
Hartley [11] considers the statistic U = W/S, where W has a distribution 
involving a scale parameter 9, S is independent of W, and s2/0% has a chi- 
square distribution with v degrees of freedom. In our problem, W = max x. 
and S = % 5s. We have taken o = 1 without loss of generality. Thus 
yY=v7vU. Let H denote the distribution function of U. It has been shown 


by Hartley that 


(3.1) u $4 (u) = vEH yo(u) - Ho(u)]. 


Hartley's method of solving the differential-difference equation (3.1) can 
be described briefly as follows. 

Regarding the degrees of freedom v as a continuous second variable taking 
on positive values, we seek a function H(v,u) of two variables for which 


(3.2) u 2Hvs4) _ yru(ve2,u) - H(v,u)] = 0 


Qu 
and 


(3.3) lim H(v,u) = F(u). 


V0 
Equation (3.2) is now converted into a partial differential equation by 
expanding the finite difference in a Taylor series. Taking the expansion 


with two more additional terms than Hartley did, we get 


4 oH Ry - 9 
5° a5 
av av 
where the arguments of H are u and v, and those of H are u and some mean 
value between v and yv+2. By introducing as new variables y = log u and 


x = log v, the partial differential equation (3.4) is transformed into 


(3.5) ae - 2 t = ze a) + 3 oh 3 “y +2 #) + 
geome 6 hen Ah. 6 as 
fp EE 10 As a5 SH 50 ho 2 Be 
iz Se 15 Shs a5 2H 205 7H 


2- < 
274 24 120 at), 


where the arguments of H are y and x, and those of H are y and some 
mean value between x and Jog(e* + 2). 


Hartley obtains an approximate solution to (3.5) for large and 


moderate values of x with the help of an iterative process. When x + ~, 

the equation (3.5) tends to the limiting form a - 2 a = 0 to which 
Hartley obtains the solution Hg(x.y) = F(eyt2%) This provides the first 
approximation. Expressed in terms of u and v, Ho depends on vv u. This 
means that, for large v, G (u) = 4) is equal to F(u) to a first 
approximation (which is independent of v). Successive closer approximations 
are obtained by Hartley iteratively. His method of solving the non-homogeneous 
partial differential equation a 2 za = y(y,x) uses the theory of 
characteristics. For details of the steps from this point on, we refer to 
Hartley's paper. We will just give here the results of Hartley for the 
three correction terms involving ae oo and Gee; and add a fourth term 
involving oe which is derived by carrying out one more iteration. The 
details of derivation are not given here. 


Let | 


m 
(3.6) G mt) = Fy) 4+ Lao 
1= 


where L; is the ith correction tern. Also, let FI) (y) denote the jth 
order derivative of Fly) w.r.t. y. Then the first four correction terms 
can be written as 
tyty) = 2 fal?) _ 03, 
1 2 1 
Loly) = az Laa'4) 1993) + 9a?) 24109, 


(3.7) L(y) = 23 (068) 7465) 4 174 (4)_47,3) , 60'2) 7, 


L(y) = 1, (15.8) - 18007) + 8300 (6) - 
360v 
184805) + 20154) ~ 90003) 4 29,2) + 48, (1) 7, 
where 


ad) = gltyy, j = 1,2,...,8, 
2 


My) - y FY), 
y) = yPFl2iy) 4 y Fy), 
oly) - VFDViy) 4 ay?Fl2)(y) 4 y FY), 
oy) = yi) (y) 4 6y3e(3)(y) 4 1y*F(2)(y) 4 y Fey), 
(3.8) oF) (y) = SEC (y) + royal (yy 4 25y36(3) (yy ¢ 
wy F(Z (yy ey (yy, 
Vy) = YE) y) + 15y5F(5) (yy + e5y4e( yy ¢ 
90y3F3)(y) 4 a1yF(2) (yy 4 y Fy), 
oy) = yy) 4 2iySr(6) cy) 4 vaoy®F 5) (y) + 
350y%F'9)(y) + sorye((y) + 6ay2F(2 (yy 4 y RLM yy, 


(8) (y) = y8r (8) (yy + 2ay7F\7) (yy + 266y%r(6)(y) 4 
(3.8 1osoy?F (5) (y) + r701yte (4) (yy + 
o66y-F (>) (y) + i27y*F (2) (y) + y FN (yy, 


cont.) 


The correction term L(y) given without the details of derivation by 
Hartley in terms of the yd) is found to be in error; stated in terms 
of the a4), his coefficient of a?) is -22/3 instead of -17. 


Now, it is easy to see that 


(Heyy = Lp gkp vets. 
(3.9) ICY) = ayy J 9 CE EEM (thet) dt, 


where H(t) is the Hermite polynomial in t of degree j. Thus the 
evaluation of Go mY) in (3.6) involves the evaluation of the integral 
in (3.9) for j = 0,1,...,8. 

Before we discuss the details of our numerical evaluation of the 
percentage points of the distribution of Y by using the approximation 
G) mY) to its distribution function G (y). we will give a brief account 


of the existing tables. 


4. EXISTING TABLES: BRIEF ACCOUNT 

As we noted earlier, several authors have discussed the distribution 
of Y and have provided relevant tables. We give a brief description of 
these tables. 

Tables in the case of v = © (i.e. the percentiles of W) have been 
given by Bechhofer [1], Gupta [6], Gupta, Nagel and Panchapakesan [8], 
and Milton [14]. Bechhofer's table gives several percentage points (ranging 
from 5% to 99.95%) to four decimal places for k = 1(1)9 and p = 0.5. 
Gupta [6] has tabulated the 75, 90, 95, 97.5, and 99 percentage points to 


three decimal places for o = 0.5 and k = 1(1)50. Gupta [6] has also given 


the values of Fi (ws 0) to five decimal places for k = 1(1)12, w = -3.5(0.1)3.5 
and seventeen selected values of p ranging from 0.1 to 0.9. The values of 

Fy (ws po) have also been provided by Milton [14] to eight decimal places for 

p = 0.5, k = 2(1)9(5)24, and w = 0.00(0.05)5.15. Finally, the table of 

Gupta, Nagel and Panchapakesan [8] gives the 75, 90, 95, 97.5 and 99 percentage 
yoints for k = 1(1)10(2)50 and the seventeen values of op chosen by Gupta 

6]. 

Percentage points of Y have been earlier tabulated by several authors. 
Pillai and Ramachandran [15] have tabulated the 95 percent points to two 
decimal places for o = 0, k = 1(1)8, and v = 3(1)10,12,14,15,16,18,20,24, 30, 
40,60,120 and ~. 

Dunnett and Sobel [4] have considered only the bivariate t (k=2) and 
tabulated the percentage points (chree decimal places) as well as the 
values of Gy) (five decimal places) for ¢ = + 0.5 and 
v = 1(1)30(3)60(15)120,150,300,600,-. Their table of the values of 


G, (y) covers y = 0.00(0.25)2.50(0.50)10.00 and the percentage levels 


k,v 
chosen for the other table are 50, 75, 90, 95 and 99. 

Dunnett [3] and Gupta and Sobel [9] have considered the case of 
o = 0.5 and k > 2. Dunnett's table gives the 95 and 99 percent points 
to two decimal places for k = 1{1)9 and » = 5(5)25 whereas in the 
Gupta-Sobel table (correct to two decimal places) k = 2,5,10(1)16,18,20(5) 
40,50; v = 15(1)20,24,30,36,40,48,60,80,100,120,360,~ and the percentage 
Tevels are 75, 90, 95, 97.5, 99, 99.75, 99.9, and 99.95. Dunnett's method 
is to numerically evaluate Gy) for three successive values of y 


differing by 0.1 such that the desired value of the probability level is 


bracketed and then determine the percentage point by inverse interpolation. 


The Gupta-Sobel table is for the percentage points of V2 Y and they use che 
Cornish-Fisher expansion with four adjustment terms. 

Krishnaiah and Armitage [13] have tabulated the 95 and 99 percent 
points for k = 1(1)10, v = 5(1)35, and o = 0(0.1)0.9. They first evaluated 
Gy (ys o) in (2.3) for y = 0.1(0.20)6.1 for the different values of p, v, 
and k using the 40-point Gauss-Hermite quadrature formula for the inner 
integral and 48-point Gauss quadrature formula for the outer integral. 


The required percentage points were then obtained by using cubic interpolation. 


5. CONSTRUCTION OF THE NEW TABLE 
The upper 100a percent point of the distribution of Y is denoted 


by y* = y*la,ky,v,o). Its numerical value is obtained as the solution of 
(5.1) Gy ny) = ]-a. 


The secant method was used to finda value of y such that 


(Sed) IG. (y)-l +a] <€ 


v,m 


with e = 1078. 


This value of y satisfying (5.2) is our solution y*. 

Now, the numerical evaluation of Go mly) involves the evaluation of 
the integrals FI) (y) given by (3.9) for j = 0,1,...,8. This was accomplished 
by using one of two different methods, labeled Method 1 and Method 2, 
depending on the value of p; the only exception is when v = © in which case 
Method 2 was used for all the selected values of p. Method 7 evaluates 


FJ) (yy using the 60-point Gauss-Hermite quadrature formula. We will denoce 


this numerically evaluated quantity by ES) (yy, To describe Method 2, letting 


(5.3) Ae ee) = 8 See 4 


and changing the variable by the transformation u = At + By, we can write 


(3.9) as 


(5.4) FD Gy) = f eeu CEM (GEN RY Fd 


The value of this integral is approximated by 


~fs 9A+4B : 
(5.5) Ne. (uy (SP o(4 BY (SF du. 


and the integration was carried out by Gauss's method over intervals of 
length D = 0.5 starting from -9 until 9A + 4B was included. Note that, for 
convenience, the approximation for FI) yy is denoted by FG) (yy for either 
method. Whatever the method, FD (yy is evaluated as a sum of the form 
pee tts y). While evaluating the terms of this sum, we used the IMSL 


subroutine for the values of the standard normal cdf. If the exponent in 


2 
-t°/2 : 
e * was less than -50 Li.e., It, | > 10), the value of FO) y) was 


taken to be zero because for the ranges of » and v in our table, the value 


10. Also, if the value of the standard 


of the termwas not more than 10° 
normal cdf involved in the expression was less than 107°, then for k » 8, 
the value of FI) (t y) was taken to be zero since the value was not more 


than 10720 for all the » values of the table. 


Now, Method 2 is more precise than Method 1; however, it is also 
more expensive, the cost ratio being 6 to 1 in some extreme cases. In | 
order to decide on the method and the number of correction terms for the 
range of selected values of o and v, some preliminary comparative studies 
were made. 

In the case of k = 19 and a = 0.01, we computed the percentage point 
y* using both methods with m = 4 correction terms for v = 15 and 120. 


The results were as follows. 


Se 


1] 


Table 1. Upper 1% point of the distribution of Y 


k=19, m= 4 

0 v Method 1 Method 2 

0.2 15 3.96044429 3.96044429 

120 3.34250119 3.34250119 

0.3 15 3.90335614 3.90335614 

20 3. 32690446 3. 32690446 

0.5 15 3.78170980 3.7817098) 

120 3. 26561715 3. 26561715 

0.6 15 3.70147425 3.70147458 

120 3.21375097 3.21375094 

0.7 15 3.59638560 3.61431834 
120 3.14096412 3.14096740 


Based on this comparison, Method 1 was used for p = 0.2(0.1)0.6, and 
Method 2 for p = 0.7(0.1)0.9. In the same case (k = 19, « = 0.01), the 
following table gives the percentage point y* for o = 0.1 using both methods 


with 4 correction terms. 


Table 2. Upper 1% point of the distribution of Y 
k = 19, m= 4, 9 = 0.1 


v Method 1 Method 2 
15 4.01116552 4.01119416 
30 3.60648913 3.60658967 
48 3.47577569 3.47577575 
60 3. 43358259 3.43358261 

120 3.35129394 3. 35129394 


Based on this, it was decided to use, in the case of o = 0.1,Method 2 


for v < 24 and Method 1 for v > 30. 


To illustrate the difference between using three and four correction 


terms, we give below the upper 1] percent point in the case of k = 17 and 


0 = 0.7 for v = 48, 60, and 120 using Method 2. 


— oa EN So ge tae oe ae 


Table 3. Upper 1% point of the distribution of Y 
k = 17, 9 = 0.7 


v 3 terms 4 terms 
48 2.94535892 2.94537526 
60 2.92086206 2.92086933 

120 2.87288348 2.87288399 


For a high value of p, three terms seem to be quite adequate for v > 60. 
For small values of p, we used three terms for v = 120 and four terms for 


other selected values of p which are all less than or equal to 60. 


We summarize below the ranges of p, k, v, anda for which the percentage 
points of Y are given in Table 4 along with information regarding the method 
(i.e.,Method lor 2) and the number of correction terms (m) used. In Table 4, 
it should be noted that P* = I-x. The entries in the table are accurate to five 


decimal places. 


The coverage of Table 4 
(i) 9 = 0.13; k = 1(1)9(2)19; a = 1-P* = 0.01,0.05,0.10,0.25 


v Method m 
15(1)20,24 2 4 
30 , 36 ,48,60 1 4 
120, T 3 


(ii) oo = 0.2(0.1)0.6; k = 1(1)9(2)19; a = 1-P* = 0.01,0.05,0.10,0.25; Method 1 


) m 

15(1) 20,24, 30, 36,48,60 4 

120, 3 
(iii) » = 0.7(0.1)0.9; k = 1(1)9(2)15; a = 1-P* = 0.05,0.10; Method 2 

v m 

15,17,20,24,36 4 

60,120,« 3 
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The entry in each case is the value of y* satisfying Gy by* oe) = Pe = |x, 


by (2.3). 


is given 


where G. (ye) 
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6. APPLICATIONS 
In this section we discuss briefly some specific applications il]lustra- 
ting the use of the tables. These applications relate to selection procedures, 


multiple comparisons, prediction intervals and tests of hypotheses. 


6.1 Selection procedures 


Let n oy be k independent normal populations with unknown 


pre 
MEANS i) 5+-+ soups respectively, and common unknown variance se Gupta 


[5], [7] proposed a procedure Ry for selecting a nonempty subset from the 
k populations so that the population associated with the largest Wy is 
included in the selected subset (this event is called a correct selection 
(CS)) with a probability at least equal to P* (1/k < P* < 1). Let 
i = 1,...,k, be the means of samples of size n from these populations. 


The rule Ry is: Select n if and only if Y; > max Y. - ds | where 5¢ is the 
1<j<k J on 


usual pooled estimator of a® based on all the k samples such that usryee 
has a yo distribution with v = k(n-1) degrees of freedom. The infimum of 
PCS, the probability of a CS, occurs when My Feet uy for all values of a > 0: 
Thus the value of d satisfying the probability requirement is given by 
Gy yld/ 25 1/2) = P*. Therefore, d = 72 y*(1-P*,k-1,k(n-1),0.5). 
Suppose the above k populations represent experimental categories and 
they are compared with a control population 9 which is N(uugs0") where ug 
is unknown. Let Yo = Xo be the mean of a sample of size m from TQ: The 
pooled estimator 5° is now based on v = k(n-1)+(m-1) degrees of freedom. 
Gupta and Sobel [10] proposed a rule Ro for selecting a subset of the k 
experimental populations such that all populations better than the control 
(i.e. those for which uy 2 ug) are included in the selected subset with a 


probability at least equal to P*. We propose a modified procedure Ro 


Select 7, if and only if XK, > X - Ds /* + u , where D > 0 is to be detennined 
to satisfy the probability requirement. It is easy to show that the infimum 


of PCS is equal to Gy, (Ds o) where p = n/(ntm). This gives D = y*(1-P*,k,v,0). 


6.2 Multiple comparisons between k treatment means and a control mean 


The problem of the simultaneous one-sided comparisons of the means of 
k treatments with that of a control group has been discussed by Dunnett 


[3]. Let X,,X,,...,X, be independent sample means based on m,n,...,n 


k 
independent observations from normal populations with unknown means 
lig ety rece stys respectively, and a common unknown variance of, Let 5¢ 
be the usual pooled unbiased estimator of of based on v = k({n-1)+(m-1) degrees of 
freedom such that vs¢/o* has a yo-distribution with v degrees of freedom. Then 
a set of k exact one-sided (lower) 100y% simultaneous confidence intervals 


for the differences Lig» i= 1,...,k, is given by 


(X;-Xp) - y*(1-y,k,v.o) S/= + 2, 


where p = n/(ntm). 


6.3 Prediction intervals to contain all of k future means 
Let YyoreeoYy be the n observations of a current sample from a normal 
population with mean u and variance 0, both unknown. Let Y be the sample 
mean and 5? be the usual unbiased estimator of Pe on v = n-l degrees of 
freedom. Let Ky r eee ok be the (unknown) sample means of k additional future 
independent samples each of size m. The x. are independent of Y and s*. An 


exact lower simultaneous 100y% prediction interval for the sample means 


Ky vee eoXy is given by ¥ - y*(1-y,k,v,0)s [i+ - where o = m/(mtn). 


; 
' 
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6.4 Test procedures for multiple comparisons 

Krishnaiah [12] proposed test procedures for the multiple comparisons 
of means (mean vectors) against one-sided alternatives under a general 
analysis of variance (multivariate analysis of variance) model. Let us 
consider k normal populations with unknown means Hy oe ee ably and a common 
unknown variance ao Consider the problem of testing 
H;: Chu = 0 simultaneously against the respective alternatives 
Ke: Cou > 0, i = 1,...,q (q < k), where x’ = Ciyoeseoty)s 
and the c; are kx] vectors of known elements such that cic, =], 
c. 1 = 0 for i = 1,...,q, and cic, =, (+ 0) for i # j = 1,....q. Here 1 
denotes a column vector whose elements are unity. It is required that 


q 
priany H. is rejected; ni 
i=] 


Hi} = ae 
Let ., i = 1,...,k, the sample means based on n independent observations 
from each population and let 52 denote the usual pooled unbiased estimator 
of of on . = k(n-1) degrees of freedom. Then, as a special case of the 
procedure proposed by Krishnaiah [12], we use the procedure: 

Yn ciZ 


j i i = = 
Reject H, (against K;) if t; “as t 


where Z = (¥).---0¥)) and to satisfying the error restriction is given 


by Ee = y*(l-a,k,vuyn). 


3] 
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